require(phyloseq)
require(tidyverse)
require(reshape2)
require(dplyr)
require(ggplot2)
require(vegan)
Load data order, factors, and create a mode (chemical, hand, non-treated) column.
ps_dmn <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/DMN_ests_16S.Rdata")
sample_data(ps_dmn)$Herbicide <- factor(sample_data(ps_dmn)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_dmn)$herb_time<-paste(sample_data(ps_dmn)$Herbicide, sample_data(ps_dmn)$Time, sep = "_")
#regroup all chemical treatments together and rerun betadiv calcs within group.
sample_data(ps_dmn)$Mode<-sample_data(ps_dmn)$Herbicide
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")
sample_data(ps_dmn)$Mode<- as.factor(values[match(sample_data(ps_dmn)$Mode, index)])
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")
sample_data(ps_dmn)$Herbicide <- as.factor(values[match(sample_data(ps_dmn)$Herbicide, index)])
ps_rare <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/HerbPt1_rare_16S.Rdata")
sample_data(ps_rare)$Herbicide <- factor(sample_data(ps_rare)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_rare)$herb_time<-paste(sample_data(ps_rare)$Herbicide, sample_data(ps_rare)$Time, sep = "_")
#regroup all chemical treatments together and rerun betadiv calcs within group.
sample_data(ps_rare)$Mode<-sample_data(ps_rare)$Herbicide
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")
sample_data(ps_rare)$Mode<- as.factor(values[match(sample_data(ps_rare)$Mode, index)])
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")
sample_data(ps_rare)$Herbicide <- as.factor(values[match(sample_data(ps_rare)$Herbicide, index)])
ps_trans <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/HerbPt1_hel_trans_16S.Rdata")
sample_data(ps_trans)$Herbicide <- factor(sample_data(ps_trans)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_trans)$herb_time<-paste(sample_data(ps_trans)$Herbicide, sample_data(ps_trans)$Time, sep = "_")
#regroup all chemical treatments together and rerun betadiv calcs within group.
sample_data(ps_trans)$Mode<-sample_data(ps_trans)$Herbicide
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")
sample_data(ps_trans)$Mode<- as.factor(values[match(sample_data(ps_trans)$Mode, index)])
index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")
sample_data(ps_trans)$Herbicide <- as.factor(values[match(sample_data(ps_trans)$Herbicide, index)])
create alpha diversity tables
alpha_div <- estimate_richness(physeq = ps_rare, measures = c("Observed", "Shannon", "Chao1"))
#pull out metadata and concatonate with alpha diversity metrics
md<-data.frame(sample_data(ps_rare))
alpha_div_md <- rownames_to_column(alpha_div, "Barcode_ID_G") %>% full_join(md)
Joining, by = "Barcode_ID_G"
alpha_div_md$Herbicide <- factor(alpha_div_md$Herbicide, levels = c("Non-Treated", "Handweeded", "Atrazine-Mesotrione", "Dicamba", "Glyphosate"))
Shannon Div plots - no significant differences among herbicide treatments at any of the three time points
ggplot(data = alpha_div_md, aes(Herbicide, Shannon, color= Herbicide)) + facet_grid(. ~ Time) + geom_boxplot() + theme_classic() + theme(axis.text.x = element_text(angle = 45, hjust = 1) )
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_Shannon.pdf")
Saving 7.29 x 4.51 in image
aov_t1<-aov(Chao1 ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T1",])
plot(aov_t1$residuals)
summary(aov_t1)
Df Sum Sq Mean Sq F value Pr(>F)
Herbicide 4 66537 16634 0.099 0.982
Residuals 51 8553352 167713
aov_t2<-aov(Chao1 ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T2",])
plot(aov_t2$residuals)
summary(aov_t2)
Df Sum Sq Mean Sq F value Pr(>F)
Herbicide 4 652901 163225 0.779 0.545
Residuals 49 10272922 209651
aov_t3<-aov(Chao1 ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T3",])
plot(aov_t3$residuals)
summary(aov_t3)
Df Sum Sq Mean Sq F value Pr(>F)
Herbicide 4 433435 108359 0.465 0.761
Residuals 50 11641251 232825
remove outliers and rare reads with less than 2 total reads
ps_dmn <- subset_samples(ps_dmn, sample_names(ps_rare) != "G166SG")
ps_rare <- subset_samples(ps_rare, sample_names(ps_rare) != "G166SG")
ps_rare_sub<-prune_taxa(taxa_sums(ps_rare) > 2, ps_rare)
ps_trans_sub<-prune_taxa(taxa_sums(ps_trans) > 0.05, ps_trans)
ordinations and adonis testing with three separate objects (i.e., dmn, rarefied, transformed). Rare taxa are removed from rarefied and transfomred to sucessfully ordinate. At this point, the transformed data will not ordinate.
ord_dmn<-ordinate(physeq = ps_dmn, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit=1000)
Run 0 stress 0.1117941
Run 1 stress 0.1117967
... Procrustes: rmse 0.00125091 max resid 0.01384803
Run 2 stress 0.1134329
Run 3 stress 0.1117728
... New best solution
... Procrustes: rmse 0.01276998 max resid 0.1385506
Run 4 stress 0.1118649
... Procrustes: rmse 0.01087397 max resid 0.1023859
Run 5 stress 0.1117959
... Procrustes: rmse 0.01364467 max resid 0.1486566
Run 6 stress 0.1115948
... New best solution
... Procrustes: rmse 0.004588589 max resid 0.05391228
Run 7 stress 0.111766
... Procrustes: rmse 0.00442201 max resid 0.05494826
Run 8 stress 0.1117957
... Procrustes: rmse 0.01142902 max resid 0.1397853
Run 9 stress 0.1115974
... Procrustes: rmse 0.0004616654 max resid 0.005151518
... Similar to previous best
Run 10 stress 0.1172805
Run 11 stress 0.1117979
... Procrustes: rmse 0.01131905 max resid 0.1393504
Run 12 stress 0.1117954
... Procrustes: rmse 0.01201186 max resid 0.1419067
Run 13 stress 0.1118117
... Procrustes: rmse 0.01155646 max resid 0.1407441
Run 14 stress 0.1117657
... Procrustes: rmse 0.004425952 max resid 0.05523005
Run 15 stress 0.1134278
Run 16 stress 0.1117966
... Procrustes: rmse 0.01230531 max resid 0.1482626
Run 17 stress 0.1115972
... Procrustes: rmse 0.0004124194 max resid 0.004505644
... Similar to previous best
Run 18 stress 0.1115948
... New best solution
... Procrustes: rmse 0.0001752304 max resid 0.00151237
... Similar to previous best
Run 19 stress 0.1117956
... Procrustes: rmse 0.01141319 max resid 0.1396299
Run 20 stress 0.1117961
... Procrustes: rmse 0.01230562 max resid 0.1477826
*** Best solution repeated 1 times
ord_rare<-ordinate(physeq = ps_rare_sub, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit=1000)
Square root transformation
Wisconsin double standardization
Run 0 stress 0.1959743
Run 1 stress 0.1958498
... New best solution
... Procrustes: rmse 0.00736199 max resid 0.08432568
Run 2 stress 0.1967565
Run 3 stress 0.1969623
Run 4 stress 0.195836
... New best solution
... Procrustes: rmse 0.004183262 max resid 0.04410732
Run 5 stress 0.1972166
Run 6 stress 0.1972004
Run 7 stress 0.1967507
Run 8 stress 0.1972118
Run 9 stress 0.1955227
... New best solution
... Procrustes: rmse 0.01439727 max resid 0.1055675
Run 10 stress 0.1967629
Run 11 stress 0.1955951
... Procrustes: rmse 0.01280822 max resid 0.1482606
Run 12 stress 0.1970466
Run 13 stress 0.1967622
Run 14 stress 0.1955933
... Procrustes: rmse 0.01278164 max resid 0.1478142
Run 15 stress 0.1958344
... Procrustes: rmse 0.01355004 max resid 0.1053769
Run 16 stress 0.1960805
Run 17 stress 0.1965852
Run 18 stress 0.1955274
... Procrustes: rmse 0.001218077 max resid 0.009525113
... Similar to previous best
Run 19 stress 0.1955259
... Procrustes: rmse 0.001157442 max resid 0.01064565
Run 20 stress 0.1955251
... Procrustes: rmse 0.0006025696 max resid 0.004705696
... Similar to previous best
*** Best solution repeated 2 times
full_ord_rare<-ggordiplots::gg_ordiplot(ord = ord_rare, groups = data.frame(sample_data(ps_rare))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
full_ord_rare$plot + theme_classic()
ord_transformed<-ordinate(physeq = ps_trans_sub, method = "NMDS", distance = "bray", trymax= 300, maxit=1000)
Run 0 stress 0.07436383
Run 1 stress 0.07356844
... New best solution
... Procrustes: rmse 0.05257797 max resid 0.5365183
Run 2 stress 0.07356423
... New best solution
... Procrustes: rmse 0.05700177 max resid 0.5049506
Run 3 stress 0.07475784
Run 4 stress 0.08233187
Run 5 stress 0.07371847
... Procrustes: rmse 0.03580052 max resid 0.2670563
Run 6 stress 0.0748452
Run 7 stress 0.07432888
Run 8 stress 0.07356475
... Procrustes: rmse 0.00101334 max resid 0.01257088
Run 9 stress 0.07351843
... New best solution
... Procrustes: rmse 0.03301646 max resid 0.277146
Run 10 stress 0.07481794
Run 11 stress 0.08269937
Run 12 stress 0.07412137
Run 13 stress 0.07564989
Run 14 stress 0.08308348
Run 15 stress 0.07434076
Run 16 stress 0.07435862
Run 17 stress 0.07790567
Run 18 stress 0.07413747
Run 19 stress 0.07388433
... Procrustes: rmse 0.05802289 max resid 0.4890175
Run 20 stress 0.0721808
... New best solution
... Procrustes: rmse 0.04835115 max resid 0.5098029
Run 21 stress 0.07432931
Run 22 stress 0.07439095
Run 23 stress 0.0741388
Run 24 stress 0.07215009
... New best solution
... Procrustes: rmse 0.001801121 max resid 0.02249514
Run 25 stress 0.07401767
Run 26 stress 0.07207119
... New best solution
... Procrustes: rmse 0.02931455 max resid 0.2673238
Run 27 stress 0.07804942
Run 28 stress 0.07293172
Run 29 stress 0.07458269
Run 30 stress 0.07630386
Run 31 stress 0.08146302
Run 32 stress 0.08008002
Run 33 stress 0.07439049
Run 34 stress 0.08424865
Run 35 stress 0.07343365
Run 36 stress 0.07356858
Run 37 stress 0.07618202
Run 38 stress 0.07475861
Run 39 stress 0.07256438
... Procrustes: rmse 0.02950194 max resid 0.267304
Run 40 stress 0.08201316
Run 41 stress 0.0866015
Run 42 stress 0.07209557
... Procrustes: rmse 0.001905361 max resid 0.02317822
Run 43 stress 0.08455037
Run 44 stress 0.08007292
Run 45 stress 0.07433954
Run 46 stress 0.07382298
Run 47 stress 0.07444042
Run 48 stress 0.07414279
Run 49 stress 0.07353217
Run 50 stress 0.0731359
Run 51 stress 0.07398746
Run 52 stress 0.07214984
... Procrustes: rmse 0.02925501 max resid 0.2670558
Run 53 stress 0.07357666
Run 54 stress 0.4158152
Run 55 stress 0.07395362
Run 56 stress 0.08199542
Run 57 stress 0.07565232
Run 58 stress 0.07207392
... Procrustes: rmse 0.001457162 max resid 0.01630803
Run 59 stress 0.07352626
Run 60 stress 0.07438306
Run 61 stress 0.08008342
Run 62 stress 0.07356382
Run 63 stress 0.07388934
Run 64 stress 0.07398744
Run 65 stress 0.0735644
Run 66 stress 0.07374194
Run 67 stress 0.07372377
Run 68 stress 0.07356866
Run 69 stress 0.08315409
Run 70 stress 0.07346563
Run 71 stress 0.07432932
Run 72 stress 0.07419749
Run 73 stress 0.07371894
Run 74 stress 0.07475716
Run 75 stress 0.07356427
Run 76 stress 0.07389694
Run 77 stress 0.07356259
Run 78 stress 0.07398762
Run 79 stress 0.07357632
Run 80 stress 0.07451023
Run 81 stress 0.07358983
Run 82 stress 0.0761558
Run 83 stress 0.07353107
Run 84 stress 0.07436762
Run 85 stress 0.07421491
Run 86 stress 0.07659576
Run 87 stress 0.08662864
Run 88 stress 0.07565735
Run 89 stress 0.07293179
Run 90 stress 0.07478225
Run 91 stress 0.07371931
Run 92 stress 0.07411583
Run 93 stress 0.07381286
Run 94 stress 0.07214946
... Procrustes: rmse 0.02928908 max resid 0.2677357
Run 95 stress 0.07214635
... Procrustes: rmse 0.002841601 max resid 0.02883108
Run 96 stress 0.08146367
Run 97 stress 0.07209518
... Procrustes: rmse 0.00213922 max resid 0.02322161
Run 98 stress 0.07436743
Run 99 stress 0.07436394
Run 100 stress 0.07475812
Run 101 stress 0.07353144
Run 102 stress 0.07347476
Run 103 stress 0.08005355
Run 104 stress 0.07631888
Run 105 stress 0.08008087
Run 106 stress 0.07634158
Run 107 stress 0.07390232
Run 108 stress 0.08005734
Run 109 stress 0.07388468
Run 110 stress 0.07451177
Run 111 stress 0.0722182
... Procrustes: rmse 0.02933714 max resid 0.2670719
Run 112 stress 0.07565095
Run 113 stress 0.073713
Run 114 stress 0.07727503
Run 115 stress 0.07360437
Run 116 stress 0.07207283
... Procrustes: rmse 0.0002854145 max resid 0.003421485
... Similar to previous best
*** Best solution repeated 1 times
Adonis testing of herbicide treatments by time point
ps_adonis<-function(physeq){
otu_tab<-data.frame(phyloseq::otu_table(physeq))
md_tab<-data.frame(phyloseq::sample_data(physeq))
if(taxa_are_rows(physeq)== T){
physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
print(anova(vegan::betadisper(physeq_dist, md_tab$Herbicide)))
vegan::adonis(physeq_dist ~ Herbicide * Time + Total_Weed_Veg , data = md_tab, permutations = 1000)
}
remove one sample with no vegetation measurement.
ps_rare_sub_57<-subset_samples(ps_rare_sub, sample_names(ps_rare_sub) != "G065SG")
ps_adonis(ps_rare_sub_57)
ps_dmn_57<-subset_samples(ps_dmn, sample_names(ps_dmn) != "G065SG")
ps_adonis(ps_dmn_57)
Ordination plots DMN
ord_t1_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)
T1_dmn<-ggordiplots::gg_ordiplot(ord = ord_t1_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_dmn$plot + theme_classic() + xlim(-0.4, 0.4) + ylim(-0.4, 0.4)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T1.pdf")
ord_t2_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T2"), method = "NMDS", distance = "bray", k=3, trymax= 100)
T2_dmn<-ggordiplots::gg_ordiplot(ord = ord_t2_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T2_dmn$plot + theme_classic()+ xlim(-0.4, 0.4) + ylim(-0.4, 0.4)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T2.pdf")
ord_t3_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T3"), method = "NMDS", distance = "bray", k=3, trymax= 100)
T3_dmn<-ggordiplots::gg_ordiplot(ord = ord_t3_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T3_dmn$plot + theme_classic()+ xlim(-0.4, 0.4) + ylim(-0.4, 0.4)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T3.pdf")
Ordination plots rarefied
ord_t1_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Square root transformation
Wisconsin double standardization
Run 0 stress 0.1810491
Run 1 stress 0.181468
... Procrustes: rmse 0.05480179 max resid 0.2537838
Run 2 stress 0.1812637
... Procrustes: rmse 0.02578262 max resid 0.1092543
Run 3 stress 0.1885435
Run 4 stress 0.1823917
Run 5 stress 0.1820307
Run 6 stress 0.1873026
Run 7 stress 0.1831937
Run 8 stress 0.1877083
Run 9 stress 0.1816042
Run 10 stress 0.1821568
Run 11 stress 0.1816494
Run 12 stress 0.1815751
Run 13 stress 0.1812635
... Procrustes: rmse 0.02573893 max resid 0.109003
Run 14 stress 0.180928
... New best solution
... Procrustes: rmse 0.02695614 max resid 0.1109664
Run 15 stress 0.1812552
... Procrustes: rmse 0.01647177 max resid 0.08102514
Run 16 stress 0.1814158
... Procrustes: rmse 0.02898979 max resid 0.1057012
Run 17 stress 0.1837645
Run 18 stress 0.1833285
Run 19 stress 0.1822759
Run 20 stress 0.1812636
... Procrustes: rmse 0.01524903 max resid 0.08041434
Run 21 stress 0.1819126
Run 22 stress 0.1837641
Run 23 stress 0.1812579
... Procrustes: rmse 0.01876556 max resid 0.09463233
Run 24 stress 0.1826156
Run 25 stress 0.1831936
Run 26 stress 0.1820917
Run 27 stress 0.181525
Run 28 stress 0.180942
... Procrustes: rmse 0.01975739 max resid 0.09890697
Run 29 stress 0.1823225
Run 30 stress 0.1854467
Run 31 stress 0.1810715
... Procrustes: rmse 0.03976112 max resid 0.1470678
Run 32 stress 0.1916469
Run 33 stress 0.1815252
Run 34 stress 0.1816037
Run 35 stress 0.1814036
... Procrustes: rmse 0.03959269 max resid 0.1481975
Run 36 stress 0.1822264
Run 37 stress 0.1825052
Run 38 stress 0.1818223
Run 39 stress 0.1810379
... Procrustes: rmse 0.01367068 max resid 0.08576246
Run 40 stress 0.1837642
Run 41 stress 0.1823669
Run 42 stress 0.1815993
Run 43 stress 0.1863375
Run 44 stress 0.1810309
... Procrustes: rmse 0.01329809 max resid 0.08407502
Run 45 stress 0.1824075
Run 46 stress 0.1970325
Run 47 stress 0.1810496
... Procrustes: rmse 0.02723117 max resid 0.1117951
Run 48 stress 0.1815987
Run 49 stress 0.1828642
Run 50 stress 0.1831933
Run 51 stress 0.1814147
... Procrustes: rmse 0.02883545 max resid 0.1054407
Run 52 stress 0.1830535
Run 53 stress 0.1812633
... Procrustes: rmse 0.01538529 max resid 0.08064122
Run 54 stress 0.1812549
... Procrustes: rmse 0.01637034 max resid 0.08068132
Run 55 stress 0.1821578
Run 56 stress 0.1881487
Run 57 stress 0.1814143
... Procrustes: rmse 0.02873776 max resid 0.1051949
Run 58 stress 0.1873034
Run 59 stress 0.1810298
... Procrustes: rmse 0.01306299 max resid 0.08278488
Run 60 stress 0.1814024
... Procrustes: rmse 0.03943623 max resid 0.1473544
Run 61 stress 0.1809421
... Procrustes: rmse 0.0198368 max resid 0.0992024
Run 62 stress 0.1810715
... Procrustes: rmse 0.03980534 max resid 0.1467033
Run 63 stress 0.182157
Run 64 stress 0.1833543
Run 65 stress 0.1810316
... Procrustes: rmse 0.01158711 max resid 0.07286033
Run 66 stress 0.1814024
... Procrustes: rmse 0.03969086 max resid 0.1472131
Run 67 stress 0.1821248
Run 68 stress 0.1823823
Run 69 stress 0.1819351
Run 70 stress 0.1816261
Run 71 stress 0.182226
Run 72 stress 0.1814204
... Procrustes: rmse 0.03583938 max resid 0.1412563
Run 73 stress 0.1819335
Run 74 stress 0.1822758
Run 75 stress 0.1819161
Run 76 stress 0.1816036
Run 77 stress 0.1814157
... Procrustes: rmse 0.02745011 max resid 0.1035174
Run 78 stress 0.1810491
... Procrustes: rmse 0.02698331 max resid 0.1109183
Run 79 stress 0.1812547
... Procrustes: rmse 0.01650078 max resid 0.08064451
Run 80 stress 0.1821955
Run 81 stress 0.1813262
... Procrustes: rmse 0.02223434 max resid 0.09864611
Run 82 stress 0.1808178
... New best solution
... Procrustes: rmse 0.009380364 max resid 0.04964183
Run 83 stress 0.1873038
Run 84 stress 0.1822261
Run 85 stress 0.182276
Run 86 stress 0.1814022
Run 87 stress 0.1832011
Run 88 stress 0.1881501
Run 89 stress 0.1907288
Run 90 stress 0.1838924
Run 91 stress 0.1818229
Run 92 stress 0.1808181
... Procrustes: rmse 0.0002255994 max resid 0.001199249
... Similar to previous best
*** Best solution repeated 1 times
T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_rare_plot<-T1_rare$plot + theme_classic() + xlim(-0.4, 0.4) + ylim(-0.4, 0.4) + guides(color=guide_legend("Treatment")) + xlab("NMDS 1") + ylab("NMDS 2")
T1_rare_plot
library(cowplot)
my_legend <- get_legend(T1_rare_plot)
library(ggpubr)
as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordinationlegend.pdf")
Saving 7.29 x 4.51 in image
as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordinationlegend.eps")
Saving 7.29 x 4.51 in image
CAP ordination plots rarefied - not used
t1_dist <- distance(subset_samples(ps_rare, Time=="T1"), method="bray") #get wUnifrac and save
t1_table<-as.matrix(dist(t1_dist)) #transform wUnifrac index
ord_t1_rare_cap <- capscale(t1_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T1"))))
T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T1_cap.pdf")
t2_dist <- distance(subset_samples(ps_rare, Time=="T2"), method="bray") #get wUnifrac and save
t2_table<-as.matrix(dist(t2_dist)) #transform wUnifrac index
ord_t2_rare_cap <- capscale(t2_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T2"))))
T2_rare<-ggordiplots::gg_ordiplot(ord = ord_t2_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T2_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T2_cap.pdf")
#G166SG identified as outlier based on plots with it included. Removed to create plot.
ps_rare <- subset_samples(ps_rare, sample_names(ps_rare) != "G166SG")
t3_dist <- distance(subset_samples(ps_rare, Time=="T3"), method="bray") #get wUnifrac and save
t3_table<-as.matrix(dist(t3_dist)) #transform wUnifrac index
ord_t3_rare_cap <- capscale(t3_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T3"))))
T3_rare<-ggordiplots::gg_ordiplot(ord = ord_t3_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T3_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T3_cap.pdf")
ggplot_build(T3_rare$plot)$data
Pairwise adonis testing
ps_pairwiseadonis<-function(physeq){
otu_tab<-data.frame(phyloseq::otu_table(physeq))
md_tab<-data.frame(phyloseq::sample_data(physeq))
if(taxa_are_rows(physeq)== T){
physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
pairwiseAdonis::pairwise.adonis(x = physeq_dist, factors = md_tab$Herbicide, p.adjust.m = "none", perm = 1000)
}
ps_t1<-subset_samples(ps_rare_sub, Time == "T1")
ps_t1<-prune_taxa(taxa_sums(ps_t1) > 2, ps_t1)
ps_t2<-subset_samples(ps_rare_sub, Time == "T2")
ps_t2<-prune_taxa(taxa_sums(ps_t2) > 2, ps_t2)
ps_t3<-subset_samples(ps_rare_sub, Time == "T3")
ps_t3<-prune_taxa(taxa_sums(ps_t3) > 2, ps_t3)
ps_pairwiseadonis(ps_t1)
ps_pairwiseadonis(ps_t2)
ps_pairwiseadonis(ps_t3)
Pairwise betadispr by treatment, time and mode
ps_betadispr<-function(physeq, groupingvar = "Groupingvar"){
otu_tab<-data.frame(phyloseq::otu_table(physeq))
md_tab<-data.frame(phyloseq::sample_data(physeq))
if(taxa_are_rows(physeq)== T){
physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
mod<-vegan::betadisper(physeq_dist, md_tab[,groupingvar])
## Perform test
print(anova(mod))
## Permutation test for F
pmod <- vegan::permutest(mod, permutations = 1000, pairwise = TRUE)
print(pmod)
print(boxplot(mod))
}
permute test of disperson
ps_betadispr(subset_samples(ps_rare_sub, Time == "T1"), groupingvar = "Mode")
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.010804 0.0054021 11.752 5.966e-05 ***
Residuals 53 0.024363 0.0004597
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.010804 0.0054021 11.752 1000 0.000999 ***
Residuals 53 0.024363 0.0004597
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Chemical Hand Non-Treated
Chemical 0.00099900 0.0030
Hand 0.00014786 0.2198
Non-Treated 0.00224710 0.23207124
$stats
[,1] [,2] [,3]
[1,] 0.3804321 0.3504290 0.3665005
[2,] 0.3977447 0.3647856 0.3788662
[3,] 0.4169569 0.3742943 0.3945778
[4,] 0.4264042 0.4009722 0.4033251
[5,] 0.4454283 0.4306733 0.4230861
$n
[1] 34 10 12
$conf
[,1] [,2] [,3]
[1,] 0.4091911 0.3562140 0.3834219
[2,] 0.4247227 0.3923745 0.4057337
$out
G046SG
0.4892527
$group
[1] 1
$names
[1] "Chemical" "Hand" "Non-Treated"
ps_betadispr(subset_samples(ps_rare_sub, Time == "T2"), groupingvar = "Mode")
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.005564 0.0027818 2.6163 0.08286 .
Residuals 51 0.054227 0.0010633
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.005564 0.0027818 2.6163 1000 0.07892 .
Residuals 51 0.054227 0.0010633
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Chemical Hand Non-Treated
Chemical 0.563437 0.0180
Hand 0.544593 0.2138
Non-Treated 0.016124 0.196261
$stats
[,1] [,2] [,3]
[1,] 0.3663187 0.3540293 0.3405369
[2,] 0.3828402 0.3788455 0.3696567
[3,] 0.3940453 0.3912101 0.3843606
[4,] 0.4182627 0.4028680 0.3904450
[5,] 0.4630042 0.4051548 0.4046331
$n
[1] 32 11 11
$conf
[,1] [,2] [,3]
[1,] 0.3841516 0.3797660 0.3744573
[2,] 0.4039391 0.4026541 0.3942639
$out
G070SG G065SG
0.5218521 0.5170988
$group
[1] 1 2
$names
[1] "Chemical" "Hand" "Non-Treated"
ps_betadispr(subset_samples(ps_rare_sub, Time == "T3"), groupingvar = "Mode")
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.000037 0.00001841 0.0189 0.9813
Residuals 51 0.049706 0.00097462
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.000037 0.00001841 0.0189 1000 0.982
Residuals 51 0.049706 0.00097462
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Chemical Hand Non-Treated
Chemical 0.95604 0.8262
Hand 0.95042 0.9421
Non-Treated 0.81436 0.93161
$stats
[,1] [,2] [,3]
[1,] 0.3660073 0.3529787 0.3641228
[2,] 0.3775720 0.3619575 0.3888203
[3,] 0.3938294 0.3833167 0.3943978
[4,] 0.4094830 0.4127647 0.4146763
[5,] 0.4548827 0.4517078 0.4523815
$n
[1] 33 11 10
$conf
[,1] [,2] [,3]
[1,] 0.3850526 0.3591128 0.3814792
[2,] 0.4026063 0.4075206 0.4073165
$out
G123SG G125SG
0.4884281 0.5089781
$group
[1] 1 2
$names
[1] "Chemical" "Hand" "Non-Treated"
ps_betadispr(subset_samples(ps_rare_sub, Mode == "Chemical"), groupingvar = "Time")
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.005636 0.00281822 3.8869 0.02381 *
Residuals 96 0.069606 0.00072506
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.005636 0.00281822 3.8869 1000 0.01798 *
Residuals 96 0.069606 0.00072506
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
T1 T2 T3
T1 0.2007992 0.0030
T2 0.1962534 0.2238
T3 0.0023253 0.1971546
$stats
[,1] [,2] [,3]
[1,] 0.3804400 0.3663109 0.3660212
[2,] 0.3977547 0.3828261 0.3776147
[3,] 0.4169648 0.3940558 0.3939165
[4,] 0.4264004 0.4182577 0.4094472
[5,] 0.4454327 0.4630282 0.4548610
$n
[1] 34 32 33
$conf
[,1] [,2] [,3]
[1,] 0.4092027 0.3841594 0.3851612
[2,] 0.4247268 0.4039521 0.4026718
$out
G046SG G070SG G123SG
0.4892518 0.5218707 0.4883976
$group
[1] 1 2 3
$names
[1] "T1" "T2" "T3"
ps_betadispr(subset_samples(ps_rare_sub, Mode == "Non-Treated"), groupingvar = "Time")
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.0019364 0.00096822 2.4191 0.1062
Residuals 30 0.0120071 0.00040024
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.0019364 0.00096822 2.4191 1000 0.1049
Residuals 30 0.0120071 0.00040024
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
T1 T2 T3
T1 0.088911 0.5215
T2 0.099548 0.0579
T3 0.496758 0.059679
$stats
[,1] [,2] [,3]
[1,] 0.3665168 0.3405704 0.3641416
[2,] 0.3788611 0.3696647 0.3888093
[3,] 0.3945750 0.3843826 0.3943811
[4,] 0.4033384 0.3904236 0.4146615
[5,] 0.4230848 0.4046057 0.4524094
$n
[1] 12 11 10
$conf
[,1] [,2] [,3]
[1,] 0.3834107 0.3744933 0.3814644
[2,] 0.4057392 0.3942719 0.4072979
$out
numeric(0)
$group
numeric(0)
$names
[1] "T1" "T2" "T3"
ps_betadispr(subset_samples(ps_rare_sub, Mode == "Hand"), groupingvar = "Time")
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.001904 0.00095187 0.5916 0.56
Residuals 29 0.046661 0.00160901
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.001904 0.00095187 0.5916 1000 0.5734
Residuals 29 0.046661 0.00160901
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
T1 T2 T3
T1 0.31868 0.3896
T2 0.27871 0.9451
T3 0.35755 0.94959
$stats
[,1] [,2] [,3]
[1,] 0.3504614 0.3540454 0.3530162
[2,] 0.3647810 0.3788398 0.3619277
[3,] 0.3742906 0.3912157 0.3833179
[4,] 0.4009630 0.4028635 0.4127718
[5,] 0.4306622 0.4051523 0.4516826
$n
[1] 10 11 11
$conf
[,1] [,2] [,3]
[1,] 0.3562127 0.3797711 0.3590964
[2,] 0.3923686 0.4026603 0.4075394
$out
G065SG G125SG
0.5171023 0.5089702
$group
[1] 2 3
$names
[1] "T1" "T2" "T3"
ps_betadispr(subset_samples(ps_rare_sub, Time == "T1"), groupingvar = "Herbicide")
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 4 0.002840 0.00070997 1.074 0.3791
Residuals 51 0.033714 0.00066106
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 4 0.002840 0.00070997 1.074 1000 0.3886
Residuals 51 0.033714 0.00066106
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Atrazine-Mesotrione Dicamba Glyphosate Handweeded Non-Treated
Atrazine-Mesotrione 0.925075 0.966034 0.086913 0.3646
Dicamba 0.930200 0.978022 0.110889 0.3996
Glyphosate 0.967636 0.975250 0.188811 0.5365
Handweeded 0.079155 0.109142 0.175310 0.2298
Non-Treated 0.369597 0.406696 0.506436 0.232071
$stats
[,1] [,2] [,3] [,4] [,5]
[1,] 0.3667119 0.3551354 0.3573486 0.3504290 0.3665005
[2,] 0.3836426 0.3774236 0.3790031 0.3647856 0.3788662
[3,] 0.4069157 0.4021368 0.3897409 0.3742943 0.3945778
[4,] 0.4088518 0.4231790 0.4188086 0.4009722 0.4033251
[5,] 0.4393221 0.4426598 0.4302367 0.4306733 0.4230861
$n
[1] 11 12 11 10 12
$conf
[,1] [,2] [,3] [,4] [,5]
[1,] 0.3949063 0.3812674 0.3707780 0.3562140 0.3834219
[2,] 0.4189251 0.4230061 0.4087038 0.3923745 0.4057337
$out
G046SG
0.4804066
$group
[1] 3
$names
[1] "Atrazine-Mesotrione" "Dicamba" "Glyphosate" "Handweeded" "Non-Treated"
ps_betadispr(subset_samples(ps_rare_sub, Time == "T2"), groupingvar = "Herbicide")
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 4 0.002843 0.00071079 0.6541 0.6268
Residuals 49 0.053249 0.00108672
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 4 0.002843 0.00071079 0.6541 1000 0.6404
Residuals 49 0.053249 0.00108672
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Atrazine-Mesotrione Dicamba Glyphosate Handweeded Non-Treated
Atrazine-Mesotrione 0.742258 0.675325 0.735265 0.4426
Dicamba 0.727584 0.313686 0.971029 0.0739
Glyphosate 0.640120 0.282922 0.399600 0.6204
Handweeded 0.734970 0.963589 0.382841 0.1918
Non-Treated 0.393767 0.097876 0.590421 0.196261
$stats
[,1] [,2] [,3] [,4] [,5]
[1,] 0.3474546 0.3612943 0.3592504 0.3540293 0.3405369
[2,] 0.3649680 0.3733844 0.3708470 0.3788455 0.3696567
[3,] 0.3810085 0.3939371 0.3751855 0.3912101 0.3843606
[4,] 0.4036468 0.4158683 0.4004167 0.4028680 0.3904450
[5,] 0.4120948 0.4539607 0.4344526 0.4051548 0.4046331
$n
[1] 11 11 10 11 11
$conf
[,1] [,2] [,3] [,4] [,5]
[1,] 0.3625824 0.3736982 0.3604113 0.3797660 0.3744573
[2,] 0.3994346 0.4141759 0.3899598 0.4026541 0.3942639
$out
G070SG G065SG
0.5064339 0.5170988
$group
[1] 1 4
$names
[1] "Atrazine-Mesotrione" "Dicamba" "Glyphosate" "Handweeded" "Non-Treated"
ps_betadispr(subset_samples(ps_rare_sub, Time == "T3"), groupingvar = "Herbicide")
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 4 0.003957 0.00098924 1.0438 0.3944
Residuals 49 0.046441 0.00094777
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 4 0.003957 0.00098924 1.0438 1000 0.3806
Residuals 49 0.046441 0.00094777
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Atrazine-Mesotrione Dicamba Glyphosate Handweeded Non-Treated
Atrazine-Mesotrione 0.369630 0.710290 0.311688 0.1049
Dicamba 0.355725 0.193806 0.764236 0.5584
Glyphosate 0.696633 0.191707 0.225774 0.0330
Handweeded 0.305819 0.741563 0.222883 0.9351
Non-Treated 0.107926 0.561747 0.033811 0.931610
$stats
[,1] [,2] [,3] [,4] [,5]
[1,] 0.3541179 0.3627381 0.3520382 0.3529787 0.3641228
[2,] 0.3639605 0.3754509 0.3633215 0.3619575 0.3888203
[3,] 0.3735777 0.3869820 0.3798037 0.3833167 0.3943978
[4,] 0.3928099 0.3954294 0.3890832 0.4127647 0.4146763
[5,] 0.4082574 0.4037548 0.4035468 0.4517078 0.4523815
$n
[1] 12 11 10 11 10
$conf
[,1] [,2] [,3] [,4] [,5]
[1,] 0.3604193 0.3774645 0.3669321 0.3591128 0.3814792
[2,] 0.3867361 0.3964996 0.3926753 0.4075206 0.4073165
$out
G131SG G123SG G125SG
0.4392663 0.4725323 0.5089781
$group
[1] 1 2 4
$names
[1] "Atrazine-Mesotrione" "Dicamba" "Glyphosate" "Handweeded" "Non-Treated"
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Glyphosate"), groupingvar = "Time")
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.0028576 0.00142878 2.0791 0.1439
Residuals 28 0.0192418 0.00068721
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.0028576 0.00142878 2.0791 1000 0.1678
Residuals 28 0.0192418 0.00068721
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
T1 T2 T3
T1 0.288711 0.0709
T2 0.260682 0.4146
T3 0.071566 0.402183
$stats
[,1] [,2] [,3]
[1,] 0.3573164 0.3592595 0.3520475
[2,] 0.3789992 0.3708523 0.3633372
[3,] 0.3897552 0.3751867 0.3798028
[4,] 0.4188308 0.4004023 0.3890938
[5,] 0.4302344 0.4344627 0.4035369
$n
[1] 11 10 10
$conf
[,1] [,2] [,3]
[1,] 0.3707798 0.3604224 0.3669337
[2,] 0.4087305 0.3899511 0.3926719
$out
G046SG
0.4803861
$group
[1] 1
$names
[1] "T1" "T2" "T3"
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Atrazine-Mesotrione"), groupingvar = "Time")
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.0020678 0.00103388 1.0923 0.348
Residuals 31 0.0293417 0.00094651
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.0020678 0.00103388 1.0923 1000 0.3716
Residuals 31 0.0293417 0.00094651
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
T1 T2 T3
T1 0.658342 0.0500
T2 0.613621 0.4595
T3 0.059278 0.439855
$stats
[,1] [,2] [,3]
[1,] 0.3667106 0.3474538 0.3540744
[2,] 0.3836487 0.3649640 0.3639284
[3,] 0.4069173 0.3810218 0.3735603
[4,] 0.4088534 0.4036548 0.3928858
[5,] 0.4393202 0.4120929 0.4082313
$n
[1] 11 11 12
$conf
[,1] [,2] [,3]
[1,] 0.3949101 0.3625900 0.3603526
[2,] 0.4189245 0.3994537 0.3867680
$out
G070SG G131SG
0.5064274 0.4391984
$group
[1] 2 3
$names
[1] "T1" "T2" "T3"
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Dicamba"), groupingvar = "Time")
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.0004756 0.00023778 0.2821 0.7561
Residuals 31 0.0261290 0.00084287
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.0004756 0.00023778 0.2821 1000 0.7572
Residuals 31 0.0261290 0.00084287
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
T1 T2 T3
T1 0.82118 0.4815
T2 0.82200 0.6314
T3 0.46827 0.62758
$stats
[,1] [,2] [,3]
[1,] 0.3551246 0.3612819 0.3627463
[2,] 0.3774302 0.3733980 0.3754452
[3,] 0.4021407 0.3938902 0.3869939
[4,] 0.4231764 0.4158655 0.3954270
[5,] 0.4426504 0.4539101 0.4037696
$n
[1] 12 11 11
$conf
[,1] [,2] [,3]
[1,] 0.3812756 0.3736592 0.3774747
[2,] 0.4230059 0.4141212 0.3965130
$out
G123SG
0.4725292
$group
[1] 3
$names
[1] "T1" "T2" "T3"
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Handweeded"), groupingvar = "Time")
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.001904 0.00095187 0.5916 0.56
Residuals 29 0.046661 0.00160901
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.001904 0.00095187 0.5916 1000 0.5774
Residuals 29 0.046661 0.00160901
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
T1 T2 T3
T1 0.31668 0.3566
T2 0.27871 0.9431
T3 0.35755 0.94959
$stats
[,1] [,2] [,3]
[1,] 0.3504614 0.3540454 0.3530162
[2,] 0.3647810 0.3788398 0.3619277
[3,] 0.3742906 0.3912157 0.3833179
[4,] 0.4009630 0.4028635 0.4127718
[5,] 0.4306622 0.4051523 0.4516826
$n
[1] 10 11 11
$conf
[,1] [,2] [,3]
[1,] 0.3562127 0.3797711 0.3590964
[2,] 0.3923686 0.4026603 0.4075394
$out
G065SG G125SG
0.5171023 0.5089702
$group
[1] 2 3
$names
[1] "T1" "T2" "T3"
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Non-Treated"), groupingvar = "Time")
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.0019364 0.00096822 2.4191 0.1062
Residuals 30 0.0120071 0.00040024
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.0019364 0.00096822 2.4191 1000 0.0989 .
Residuals 30 0.0120071 0.00040024
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
T1 T2 T3
T1 0.090909 0.5145
T2 0.099548 0.0639
T3 0.496758 0.059679
$stats
[,1] [,2] [,3]
[1,] 0.3665168 0.3405704 0.3641416
[2,] 0.3788611 0.3696647 0.3888093
[3,] 0.3945750 0.3843826 0.3943811
[4,] 0.4033384 0.3904236 0.4146615
[5,] 0.4230848 0.4046057 0.4524094
$n
[1] 12 11 10
$conf
[,1] [,2] [,3]
[1,] 0.3834107 0.3744933 0.3814644
[2,] 0.4057392 0.3942719 0.4072979
$out
numeric(0)
$group
numeric(0)
$names
[1] "T1" "T2" "T3"
ps_betadispr(ps_rare_sub, groupingvar = "Herbicide")
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 4 0.003022 0.00075548 0.8548 0.4926
Residuals 159 0.140523 0.00088379
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 4 0.003022 0.00075548 0.8548 1000 0.5185
Residuals 159 0.140523 0.00088379
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Atrazine-Mesotrione Dicamba Glyphosate Handweeded Non-Treated
Atrazine-Mesotrione 0.29870 0.73926 0.46454 0.8711
Dicamba 0.27184 0.11988 0.90010 0.1299
Glyphosate 0.73022 0.10681 0.28272 0.7942
Handweeded 0.45118 0.88577 0.27400 0.3097
Non-Treated 0.87491 0.12162 0.79878 0.31930
$stats
[,1] [,2] [,3] [,4] [,5]
[1,] 0.3635263 0.3726987 0.3644354 0.3712730 0.3653358
[2,] 0.3806232 0.3917390 0.3838596 0.3907178 0.3894408
[3,] 0.3980967 0.4079200 0.3980853 0.3985758 0.4017338
[4,] 0.4199957 0.4233531 0.4072361 0.4212779 0.4126745
[5,] 0.4700021 0.4662951 0.4375529 0.4411483 0.4390542
$n
[1] 34 34 31 32 33
$conf
[,1] [,2] [,3] [,4] [,5]
[1,] 0.3874280 0.3993536 0.3914516 0.3900402 0.3953436
[2,] 0.4087653 0.4164864 0.4047190 0.4071115 0.4081241
$out
G070SG G123SG G046SG G076SG G065SG G125SG G151SG
0.5184955 0.4995460 0.4894611 0.4554918 0.5311635 0.5552741 0.4631531
$group
[1] 1 2 3 3 4 4 5
$names
[1] "Atrazine-Mesotrione" "Dicamba" "Glyphosate" "Handweeded" "Non-Treated"
ps_betadispr(ps_rare_sub, groupingvar = "Mode")
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.001682 0.00084096 0.9478 0.3897
Residuals 161 0.142847 0.00088725
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.001682 0.00084096 0.9478 1000 0.3746
Residuals 161 0.142847 0.00088725
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Chemical Hand Non-Treated
Chemical 0.97502 0.1279
Hand 0.97790 0.2827
Non-Treated 0.13961 0.31930
$stats
[,1] [,2] [,3]
[1,] 0.3678724 0.3712730 0.3653358
[2,] 0.3914997 0.3907178 0.3894408
[3,] 0.4069683 0.3985758 0.4017338
[4,] 0.4198506 0.4212779 0.4126745
[5,] 0.4591917 0.4411483 0.4390542
$n
[1] 99 32 33
$conf
[,1] [,2] [,3]
[1,] 0.4024663 0.3900402 0.3953436
[2,] 0.4114703 0.4071115 0.4081241
$out
G046SG G070SG G087SG G123SG G131SG G065SG G125SG G151SG
0.4911052 0.5250529 0.4726503 0.4968193 0.4692099 0.5311635 0.5552741 0.4631531
$group
[1] 1 1 1 1 1 2 2 3
$names
[1] "Chemical" "Hand" "Non-Treated"
ps_betadispr(ps_rare_sub, groupingvar = "Time")
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 2 0.001119 0.00055970 0.6361 0.5307
Residuals 161 0.141654 0.00087984
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 1000
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.001119 0.00055970 0.6361 1000 0.5455
Residuals 161 0.141654 0.00087984
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
T1 T2 T3
T1 0.29271 0.3367
T2 0.29177 0.9800
T3 0.31167 0.98337
$stats
[,1] [,2] [,3]
[1,] 0.3751393 0.3651628 0.3656380
[2,] 0.3983422 0.3889047 0.3843976
[3,] 0.4102601 0.3985027 0.4009712
[4,] 0.4266630 0.4186156 0.4162396
[5,] 0.4483672 0.4529592 0.4619116
$n
[1] 56 54 54
$conf
[,1] [,2] [,3]
[1,] 0.4042806 0.3921145 0.3941248
[2,] 0.4162397 0.4048909 0.4078175
$out
G046SG G065SG G070SG G087SG G123SG G125SG
0.4859161 0.5363026 0.5212376 0.4661354 0.4943452 0.5607126
$group
[1] 1 2 2 2 3 3
$names
[1] "T1" "T2" "T3"
beta_boxplot<-function (physeq, method = "bray", group)
{
require("phyloseq")
require("ggplot2")
group2samp <- list()
group_list <- get_variable(sample_data(physeq), group)
for (groups in levels(group_list)) {
target_group <- which(group_list == groups)
group2samp[[groups]] <- sample_names(physeq)[target_group]
}
beta_div_dist <- phyloseq::distance(physeq = physeq, method = method)
beta_div_dist <- as(beta_div_dist, "matrix")
dist_df <- data.frame()
counter <- 1
for (groups in names(group2samp)) {
sub_dist <- beta_div_dist[group2samp[[groups]], group2samp[[groups]]]
no_samp_col <- ncol(sub_dist)
no_samp_row <- nrow(sub_dist)
for (cols in seq(no_samp_col)) {
if (cols > 1) {
for (rows in seq((cols - 1))) {
dist_df[counter, "sample_pair"] <- paste0(colnames(sub_dist)[cols],
"-", rownames(sub_dist)[rows])
dist_df[counter, "group"] <- groups
dist_df[counter, "beta_div_method"] <- method
dist_df[counter, "beta_div_value"] <- sub_dist[rows,
cols]
counter = counter + 1
}
}
}
}
plot_boxplot <- ggplot(data = dist_df, aes(x = group, y = beta_div_value,
color = group)) + geom_boxplot(outlier.shape = NA) +
geom_jitter() + theme_bw() + xlab(group) + ylab(method) +
theme(axis.text.x = element_text(angle = 45, vjust = 1,
hjust = 1))
list_Out <- list(data = dist_df, plot = plot_boxplot)
return(list_Out)
}
box and whisker plots of pairwise distance within group distances
#remotes::install_github("antonioggsousa/micrUBIfuns")
#library(micrUBIfuns)
T1_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T1"), method = "bray", group = "Herbicide")
T1_beta_plot <- T1_beta$plot
T1_beta_plot <- T1_beta_plot + theme_classic()+ guides(color=guide_legend("Treatment")) + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75)
T1_beta_plot
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
my_legend <- get_legend(T1_beta_plot)
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_beta_legend.pdf")
Saving 7.29 x 4.51 in image
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_beta_legend.eps")
Saving 7.29 x 4.51 in image
T1_beta_plot<-T1_beta_plot+ theme(legend.position = "none")
T1_beta_plot
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T1_rare_withingroup_beta.pdf")
Saving 7.29 x 4.51 in image
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T1_rare_withingroup_beta.eps")
Saving 7.29 x 4.51 in image
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
T1_beta_df<- T1_beta$data
T1_betamod<-aov(formula = beta_div_value ~ group ,data = T1_beta_df)
summary(T1_betamod)
Df Sum Sq Mean Sq F value Pr(>F)
group 4 0.02391 0.005978 5.4 0.000334 ***
Residuals 282 0.31223 0.001107
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
TukeyHSD(x = T1_betamod, which = "group")
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ group, data = T1_beta_df)
$group
diff lwr upr p adj
Dicamba-Atrazine-Mesotrione -0.001507071 -0.018186144 0.015172002 0.9991606
Glyphosate-Atrazine-Mesotrione 0.001896970 -0.015523754 0.019317693 0.9982525
Handweeded-Atrazine-Mesotrione -0.023576431 -0.041939486 -0.005213376 0.0044684
Non-Treated-Atrazine-Mesotrione -0.013193939 -0.029873012 0.003485134 0.1934686
Glyphosate-Dicamba 0.003404040 -0.013275033 0.020083113 0.9805780
Handweeded-Dicamba -0.022069360 -0.039730381 -0.004408340 0.0061800
Non-Treated-Dicamba -0.011686869 -0.027589741 0.004216003 0.2601887
Handweeded-Glyphosate -0.025473401 -0.043836456 -0.007110346 0.0015994
Non-Treated-Glyphosate -0.015090909 -0.031769982 0.001588164 0.0971493
Non-Treated-Handweeded 0.010382492 -0.007278529 0.028043512 0.4896459
T2_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T2"), method = "bray", group = "Herbicide")
T2_beta_plot <- T2_beta$plot
T2_beta_plot <- T2_beta_plot+ theme_classic() + theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + ggtitle("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75)
T2_beta_plot
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T2_rare_withingroup_beta.pdf")
Saving 7.29 x 4.51 in image
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T2_rare_withingroup_beta.eps")
Saving 7.29 x 4.51 in image
T2_beta_df<- T2_beta$data
T2_betamod<-aov(formula = beta_div_value ~ group ,data = T2_beta_df)
summary(T2_betamod)
Df Sum Sq Mean Sq F value Pr(>F)
group 4 0.0328 0.008212 5.471 0.000303 ***
Residuals 260 0.3902 0.001501
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
TukeyHSD(x = T2_betamod, which = "group")
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ group, data = T2_beta_df)
$group
diff lwr upr p adj
Dicamba-Atrazine-Mesotrione 0.004587879 -0.015705828 0.024881585 0.9716343
Glyphosate-Atrazine-Mesotrione -0.011421549 -0.032812994 0.009969896 0.5850476
Handweeded-Atrazine-Mesotrione 0.007709091 -0.012584615 0.028002797 0.8348219
Non-Treated-Atrazine-Mesotrione -0.022036364 -0.042330070 -0.001742657 0.0257655
Glyphosate-Dicamba -0.016009428 -0.037400872 0.005382017 0.2426725
Handweeded-Dicamba 0.003121212 -0.017172494 0.023414919 0.9933115
Non-Treated-Dicamba -0.026624242 -0.046917949 -0.006330536 0.0034255
Handweeded-Glyphosate 0.019130640 -0.002260805 0.040522085 0.1039333
Non-Treated-Glyphosate -0.010614815 -0.032006260 0.010776630 0.6518171
Non-Treated-Handweeded -0.029745455 -0.050039161 -0.009451748 0.0007049
T3_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T3"), method = "bray", group = "Herbicide")
T3_beta$plot #+ scale_color_manual(values = c("#F8766D", "#A3A500", "#00BF7D", "#00B0F6", "#E76BF3")) +
T3_beta_plot <- T3_beta$plot
T3_beta_plot <- T3_beta_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + ggtitle("")
T3_beta_plot <-T3_beta_plot + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T3_rare_withingroup_beta.pdf")
Saving 7.29 x 4.51 in image
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T3_rare_withingroup_beta.eps")
Saving 7.29 x 4.51 in image
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
T3_beta_df<- T3_beta$data
T3_betamod<-aov(formula = beta_div_value ~ group ,data = T3_beta_df)
summary(T3_betamod)
Df Sum Sq Mean Sq F value Pr(>F)
group 4 0.0477 0.011924 10.39 7.93e-08 ***
Residuals 261 0.2994 0.001147
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
TukeyHSD(x = T3_betamod, which = "group")
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ group, data = T3_beta_df)
$group
diff lwr upr p adj
Dicamba-Atrazine-Mesotrione 0.0183969697 0.001410694 0.035383246 0.0263568
Glyphosate-Atrazine-Mesotrione -0.0007666667 -0.018752976 0.017219643 0.9999574
Handweeded-Atrazine-Mesotrione 0.0276454545 0.010659179 0.044631731 0.0001130
Non-Treated-Atrazine-Mesotrione 0.0315000000 0.013513690 0.049486310 0.0000250
Glyphosate-Dicamba -0.0191636364 -0.037864911 -0.000462362 0.0415656
Handweeded-Dicamba 0.0092484848 -0.008493102 0.026990072 0.6075977
Non-Treated-Dicamba 0.0131030303 -0.005598244 0.031804305 0.3068486
Handweeded-Glyphosate 0.0284121212 0.009710847 0.047113396 0.0003917
Non-Treated-Glyphosate 0.0322666667 0.012652605 0.051880729 0.0000917
Non-Treated-Handweeded 0.0038545455 -0.014846729 0.022555820 0.9798083
library(ggpubr)
ggarrange(T1_beta_plot, T2_beta_plot, T3_beta_plot, ncol = 1)
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_within_group_beta.pdf", width = 5, height = 10)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_within_group_beta.eps", width = 5, height = 10)
Examination of dissimliarity across time points by treatment and then again by all chemical treatments combined.
T1_beta_df$Time<-"T1"
T2_beta_df$Time<-"T2"
T3_beta_df$Time<-"T3"
beta_div_T1_T2_T3 <- rbind(T1_beta_df, T2_beta_df, T3_beta_df)
beta_anova<-function(data, Herbicide = "Herbicide"){
df_sub<- data %>% filter(group == Herbicide)
mod<-aov(beta_div_value ~ Time, data = df_sub)
print(summary(mod))
print(TukeyHSD(mod, "Time"))
boxplot(df_sub$beta_div_value ~ df_sub$Time)
}
beta_anova(beta_div_T1_T2_T3, Herbicide = "Non-Treated")
Df Sum Sq Mean Sq F value Pr(>F)
Time 2 0.02547 0.01274 17.94 9.08e-08 ***
Residuals 163 0.11573 0.00071
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ Time, data = df_sub)
$Time
diff lwr upr p adj
T2-T1 -0.01621212 -0.027718988 -0.004705255 0.0030431
T3-T1 0.01578788 0.003603568 0.027972190 0.0071575
T3-T2 0.03200000 0.019331357 0.044668643 0.0000000
beta_anova(beta_div_T1_T2_T3, Herbicide = "Handweeded")
Df Sum Sq Mean Sq F value Pr(>F)
Time 2 0.01713 0.008567 4.195 0.0168 *
Residuals 152 0.31041 0.002042
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ Time, data = df_sub)
$Time
diff lwr upr p adj
T2-T1 0.02391582 0.0024158024 0.04541585 0.0251783
T3-T1 0.02231582 0.0008158024 0.04381585 0.0399463
T3-T2 -0.00160000 -0.0219967123 0.01879671 0.9811772
beta_anova(beta_div_T1_T2_T3, Herbicide = "Dicamba")
Df Sum Sq Mean Sq F value Pr(>F)
Time 2 0.00273 0.001366 0.977 0.378
Residuals 173 0.24171 0.001397
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ Time, data = df_sub)
$Time
diff lwr upr p adj
T2-T1 -0.001274747 -0.01740797 0.014858477 0.9809505
T3-T1 -0.009002020 -0.02513524 0.007131204 0.3864971
T3-T2 -0.007727273 -0.02457788 0.009123330 0.5252045
beta_anova(beta_div_T1_T2_T3, Herbicide = "Atrazine-Mesotrione")
Df Sum Sq Mean Sq F value Pr(>F)
Time 2 0.02773 0.013867 11.41 2.22e-05 ***
Residuals 173 0.21026 0.001215
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ Time, data = df_sub)
$Time
diff lwr upr p adj
T2-T1 -0.007369697 -0.02308574 0.008346347 0.5100816
T3-T1 -0.028906061 -0.04395303 -0.013859094 0.0000310
T3-T2 -0.021536364 -0.03658333 -0.006489397 0.0025367
beta_anova(beta_div_T1_T2_T3, Herbicide = "Glyphosate")
Df Sum Sq Mean Sq F value Pr(>F)
Time 2 0.02597 0.012985 14.9 1.33e-06 ***
Residuals 142 0.12374 0.000871
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = beta_div_value ~ Time, data = df_sub)
$Time
diff lwr upr p adj
T2-T1 -0.02068822 -0.03474249 -0.006633939 0.0018742
T3-T1 -0.03156970 -0.04562397 -0.017515421 0.0000012
T3-T2 -0.01088148 -0.02562173 0.003858768 0.1909546
#regroup all chemical treatments together and rerun betadiv calcs within group.
sample_data(ps_rare)$Mode<-sample_data(ps_rare)$Herbicide
index <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Handweeded", "Non-Treated")
sample_data(ps_rare)$Mode<- as.factor(values[match(sample_data(ps_rare)$Mode, index)])
#+ scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3"))
T1_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T1"), method = "bray", group = "Mode")
T1_beta_chemical_combined_plot <- T1_beta_chemical_combined$plot
T1_beta_chemical_combined_plot<- T1_beta_chemical_combined_plot + theme_classic() + guides(color=guide_legend("Treatment")) + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3"))
T1_beta_chemical_combined_plot
Warning: Removed 2 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 2 rows containing missing values (`geom_point()`).
my_legend <- get_legend(T1_beta_chemical_combined_plot)
Warning: Removed 2 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 2 rows containing missing values (`geom_point()`).
as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_beta_combined_legend.pdf")
Saving 7.29 x 4.51 in image
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_beta_combined_legend.eps")
Saving 7.29 x 4.51 in image
T1_beta_chemical_combined_plot<-T1_beta_chemical_combined_plot+ theme(legend.position = "none")
T1_beta_chemical_combined_plot
Warning: Removed 2 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 2 rows containing missing values (`geom_point()`).
T2_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T2"), method = "bray", group = "Mode")
T2_beta_chemical_combined_plot <- T2_beta_chemical_combined$plot
T2_beta_chemical_combined_plot<- T2_beta_chemical_combined_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3"))
T2_beta_chemical_combined_plot
T3_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T3"), method = "bray", group = "Mode")
T3_beta_chemical_combined_plot <- T3_beta_chemical_combined$plot
T3_beta_chemical_combined_plot<- T3_beta_chemical_combined_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3"))
T3_beta_chemical_combined_plot
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
ggarrange(T1_beta_chemical_combined_plot, T2_beta_chemical_combined_plot, T3_beta_chemical_combined_plot, ncol = 1)
Warning: Removed 2 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 2 rows containing missing values (`geom_point()`).
Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).
Warning: Removed 1 rows containing missing values (`geom_point()`).
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_within_group_beta_chemical_combined.pdf", width = 5, height = 10)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_within_group_beta_chemical_combined.eps", width = 5, height = 10)
T1_beta_df_chemical_combined <- T1_beta_chemical_combined$data
T2_beta_df_chemical_combined<- T2_beta_chemical_combined$data
T3_beta_df_chemical_combined<- T3_beta_chemical_combined$data
T1_beta_df_chemical_combined$Time<-"T1"
T2_beta_df_chemical_combined$Time<-"T2"
T3_beta_df_chemical_combined$Time<-"T3"
m1<-aov(beta_div_value ~ group, data = T1_beta_df_chemical_combined)
summary(m1)
TukeyHSD(m1, "group")
boxplot(beta_div_value ~ group, data = T1_beta_df_chemical_combined)
m2<-aov(beta_div_value ~ group, data = T2_beta_df_chemical_combined)
summary(m2)
TukeyHSD(m2, "group")
boxplot(beta_div_value ~ group, data = T2_beta_df_chemical_combined)
m3<-aov(beta_div_value ~ group, data = T3_beta_df_chemical_combined)
summary(m3)
TukeyHSD(m3, "group")
boxplot(beta_div_value ~ group, data = T3_beta_df_chemical_combined)
beta_div_chemical_combined_T1_T2_T3 <- rbind(T1_beta_df_chemical_combined, T2_beta_df_chemical_combined, T3_beta_df_chemical_combined)
beta_anova(beta_div_chemical_combined_T1_T2_T3, Herbicide = "Chemical")
beta_anova(beta_div_chemical_combined_T1_T2_T3, Herbicide = "Hand")
beta_anova(beta_div_chemical_combined_T1_T2_T3, Herbicide = "Non-Treated")
treatment to control
plotDistances = function(p, m, s, d) {
# calc distances
wu = phyloseq::distance(p, m)
wu.m = melt(as.matrix(wu))
# remove self-comparisons
wu.m = wu.m %>%
filter(as.character(Var1) != as.character(Var2)) %>%
mutate_if(is.factor,as.character)
# get sample data (S4 error OK and expected)
sd = data.frame(sample_data(p)) %>%
select(s, d) %>%
mutate_if(is.factor,as.character)
sd$Herbicide <- factor(sd$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
# combined distances with sample data
colnames(sd) = c("Var1", "Type1")
wu.sd = left_join(wu.m, sd, by = "Var1")
colnames(sd) = c("Var2", "Type2")
wu.sd = left_join(wu.sd, sd, by = "Var2")
#remove this line to plot all comparisons.
#wu.sd = wu.sd %>% filter(Type1 == "Hand" | Type1 == "Non-Treated")
# plot
ggplot(wu.sd, aes(x = Type2, y = value)) +
theme_bw() +
geom_point() +
geom_boxplot(aes(color = ifelse(Type1 == Type2, "red", "black"))) +
scale_color_identity() +
facet_wrap(~ Type1, scales = "free_x") +
theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) +
ggtitle(paste0("Distance Metric = ", m))
}
a<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T1"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
a <- a + ggtitle("Time 1 Bray-Curtis Dissimlarities")
#ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T1_rare_allgroup_beta.pdf")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T1_rare_allgroup_beta_multicomparison.pdf")
b<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T2"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
b <-b + ggtitle("Time 2 Bray-Curtis Dissimlarities")
#ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T2_rare_allgroup_beta.pdf")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T2_rare_allgroup_beta_multicomparison.pdf")
c<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T3"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
c<- c + ggtitle("Time 3 Bray-Curtis Dissimlarities")
#ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T3_rare_allgroup_beta.pdf")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T3_rare_allgroup_beta_multicomparison.pdf")
library(ggpubr)
ggarrange(a, b, c, ncol = 1)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_allgroup_beta_multi_comparison.pdf", width = 7, height = 10)
Taxon abundance bar plot
#create super long color vector
col_vector <- c("#000000", "#FFFF00", "#1CE6FF", "#FF34FF", "#FF4A46", "#008941", "#006FA6", "#A30059",
"#FFDBE5", "#7A4900", "#0000A6", "#63FFAC", "#B79762", "#004D43", "#8FB0FF", "#997D87",
"#5A0007", "#809693", "#FEFFE6", "#1B4400", "#4FC601", "#3B5DFF", "#4A3B53", "#FF2F80",
"#61615A", "#BA0900", "#6B7900", "#00C2A0", "#FFAA92", "#FF90C9", "#B903AA", "#D16100",
"#DDEFFF", "#000035", "#7B4F4B", "#A1C299", "#300018", "#0AA6D8", "#013349", "#00846F",
"#372101", "#FFB500", "#C2FFED", "#A079BF", "#CC0744", "#C0B9B2", "#C2FF99", "#001E09",
"#00489C", "#6F0062", "#0CBD66", "#EEC3FF", "#456D75", "#B77B68", "#7A87A1", "#788D66",
"#885578", "#FAD09F", "#FF8A9A", "#D157A0", "#BEC459", "#456648", "#0086ED", "#886F4C",
"#34362D", "#B4A8BD", "#00A6AA", "#452C2C", "#636375", "#A3C8C9", "#FF913F", "#938A81",
"#575329", "#00FECF", "#B05B6F", "#8CD0FF", "#3B9700", "#04F757", "#C8A1A1", "#1E6E00",
"#7900D7", "#A77500", "#6367A9", "#A05837", "#6B002C", "#772600", "#D790FF", "#9B9700",
"#549E79", "#FFF69F", "#201625", "#72418F", "#BC23FF", "#99ADC0", "#3A2465", "#922329",
"#5B4534", "#FDE8DC", "#404E55", "#0089A3", "#CB7E98", "#A4E804", "#324E72", "#6A3A4C",
"#83AB58", "#001C1E", "#D1F7CE", "#004B28", "#C8D0F6", "#A3A489", "#806C66", "#222800",
"#BF5650", "#E83000", "#66796D", "#DA007C", "#FF1A59", "#8ADBB4", "#1E0200", "#5B4E51",
"#C895C5", "#320033", "#FF6832", "#66E1D3", "#CFCDAC", "#D0AC94", "#7ED379", "#012C58",
"#7A7BFF", "#D68E01", "#353339", "#78AFA1", "#FEB2C6", "#75797C", "#837393", "#943A4D",
"#B5F4FF", "#D2DCD5", "#9556BD", "#6A714A", "#001325", "#02525F", "#0AA3F7", "#E98176",
"#DBD5DD", "#5EBCD1", "#3D4F44", "#7E6405", "#02684E", "#962B75", "#8D8546", "#9695C5",
"#E773CE", "#D86A78", "#3E89BE", "#CA834E", "#518A87", "#5B113C", "#55813B", "#E704C4",
"#00005F", "#A97399", "#4B8160", "#59738A", "#FF5DA7", "#F7C9BF", "#643127", "#513A01",
"#6B94AA", "#51A058", "#A45B02", "#1D1702", "#E20027", "#E7AB63", "#4C6001", "#9C6966",
"#64547B", "#97979E", "#006A66", "#391406", "#F4D749", "#0045D2", "#006C31", "#DDB6D0",
"#7C6571", "#9FB2A4", "#00D891", "#15A08A", "#BC65E9", "#FFFFFE", "#C6DC99", "#203B3C",
"#671190", "#6B3A64", "#F5E1FF", "#FFA0F2", "#CCAA35", "#374527", "#8BB400", "#797868",
"#C6005A", "#3B000A", "#C86240", "#29607C", "#402334", "#7D5A44", "#CCB87C", "#B88183",
"#AA5199", "#B5D6C3", "#A38469", "#9F94F0", "#A74571", "#B894A6", "#71BB8C", "#00B433",
"#789EC9", "#6D80BA", "#953F00", "#5EFF03", "#E4FFFC", "#1BE177", "#BCB1E5", "#76912F",
"#003109", "#0060CD", "#D20096", "#895563", "#29201D", "#5B3213", "#A76F42", "#89412E",
"#1A3A2A", "#494B5A", "#A88C85", "#F4ABAA", "#A3F3AB", "#00C6C8", "#EA8B66", "#958A9F",
"#BDC9D2", "#9FA064", "#BE4700", "#658188", "#83A485", "#453C23", "#47675D", "#3A3F00",
"#061203", "#DFFB71", "#868E7E", "#98D058", "#6C8F7D", "#D7BFC2", "#3C3E6E", "#D83D66",
"#2F5D9B", "#6C5E46", "#D25B88", "#5B656C", "#00B57F", "#545C46", "#866097", "#365D25",
"#252F99", "#00CCFF", "#674E60", "#FC009C", "#92896B")
phylumGlommed <- tax_glom(ps_rare, "Phylum")
#t1
phylumGlommed_herb_t1 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T1"), group = "Herbicide")
phylumGlommed_herb_t1 <- transform_sample_counts(phylumGlommed_herb_t1, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t1)$Herbicide <- factor(sample_data(phylumGlommed_herb_t1)$Herbicide, levels = c(1, 2, 3, 4, 5),
labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
plot_bar(phylumGlommed_herb_t1, x = "Herbicide", fill = "Phylum") + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) +
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_Taxon_barplot_t1.pdf")
#t2
phylumGlommed_herb_t2 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T2"), group = "Herbicide")
phylumGlommed_herb_t2 <- transform_sample_counts(phylumGlommed_herb_t2, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t2)$Herbicide <- factor(sample_data(phylumGlommed_herb_t2)$Herbicide, levels = c(1, 2, 3, 4, 5),
labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
plot_bar(phylumGlommed_herb_t2, x = "Herbicide", fill = "Phylum") + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) +
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/16S_Taxon_barplot_t2.pdf")
#t3
phylumGlommed_herb_t3 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T3"), group = "Herbicide")
phylumGlommed_herb_t3 <- transform_sample_counts(phylumGlommed_herb_t3, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t3)$Herbicide <- factor(sample_data(phylumGlommed_herb_t3)$Herbicide, levels = c(1, 2, 3, 4, 5),
labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
plot_bar(phylumGlommed_herb_t3, x = "Herbicide", fill = "Phylum") + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) +
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/16S_Taxon_barplot_t3.pdf")
Combined herbicide and time bar plot for exploration
sample_data(ps_rare)$herb_time<-paste(sample_data(ps_rare)$Herbicide, sample_data(ps_rare)$Time, sep = "_")
ps_rare_for_barplot <- prune_taxa(taxa_sums(ps_rare) > 50, ps_rare)
plot_bar(ps_rare_for_barplot, x= "herb_time", fill = "Family") + scale_fill_manual(values = col_vector) + geom_bar(stat="identity")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_BarPlot_Herbicide_Time.pdf", width = 20, height = 11)
Linear modeling of abundant taxa.
sig_mods
NULL